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Abstract 

We report on the existence and phenomenology of type II critical collapse 
within the one-parameter family of SU(2) a-models coupled to gravity. Nu- 
merical investigations in spherical symmetry show discretely self-similar (DSS) 
behavior at the threshold of black hole formation for values of the dimen- 
sionless coupling constant r\ ranging from 0.2 to 100; at 0.18 we see small 
deviations from DSS. While the echoing period A of the critical solution rises 
sharply towards the lower limit of this range, the characteristic mass scaling 
has a critical exponent 7 which is almost independent of 77, asymptoting to 
0.1185 db 0.0005 at large rj. We also hnd critical scaling of the scalar cur- 
vature for near-critical initial data. Our numerical results are based on an 
outgoing-null-cone formulation of the Einstein-matter equations, specialized 
to spherical symmetry. Our numerically computed initial-data critical pa- 
rameters p* show 2nd order convergence with the grid resolution, and after 
compensating for this variation in p* , our individual evolutions are uniformly 
2nd order convergent even very close to criticality. 
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I. INTRODUCTION 



Since the numerical investigation of dynamical behavior of a massless scalar field under 
the influence of its gravitational forces by Choptuik (Ref. [1]), critical behavior has been 
observed in a number of different matter models coupled to gravity. In the context of 
type II critical collapse, these models have in common that at the threshold of black hole 
formation their dynamics show a universal characteristic approach to either a discretely 
(DSS) or continuously (CSS) self-similar solution. 

Nonlinear a-fields provide particularly interesting models to study the dynamics of grav- 
itating self-interacting matter in general relativity. Besides their applications in physics (see 
e.g. Ref. [2]), they have a simple geometrical interpretation as harmonic maps, which have 
been extensively studied in the mathematical literature (see e.g. Refs. [3,4]). 

Recently Bizoh et al. (Refs. [5,6]) and also independently Liebling et al. (Ref. [7]) 
have observed critical (threshold) behavior for non-gravitating systems: The transition be- 
tween globally regular time evolution and singularity formation for the SU(2) cr-model on 
Minkowski background. It was shown by Bizoh (Ref. [5]) that this system admits a count- 
ably infinite family of CSS solutions. The stable ground state is the endpoint of singular 
evolution for supercritical initial data, while the first excitation, which has one unstable 
mode, plays the role of the critical (CSS) solution. 

The interesting question arises of what happens if gravity is added to this system. The 
gravitating SU(2) a-model is a family of theories with a dimensionless parameter rj, which 
acts as a coupling constant (for r\ = gravity decouples from the field). It was argued in 
Ref. [6] that the singularity formation in flat space might not be relevant for black hole 
formation when gravity is active, since the CSS blowup excludes the concentration of energy 
at the singularity. Since no asymptotically flat solitonic configurations exist (Ref. [8]), this 
suggests that the only alternative to dispersion or collapse to a black hole is the formation 
of a naked singularity. Here we focus on critical phenomena at the threshold of black hole 
formation. As Bizoh et al. have pointed out (Ref. [9]), criticality is expected to depend on 
the coupling constant rj. If so, does the system show discrete or continuous self-similarity? 
And in which way do critical phenomena depend on the coupling? 

In this paper we present results from a numerical study of the dynamical evolution for 
the SU(2) nonlinear cr-model coupled to gravity in spherical symmetry. Our code uses a 
characteristic formulation, specialized to the spherical symmetry. Initial data are specified 
on an outgoing null cone with vertex at the center of symmetry. The discretized field 
equations are used to evolve the matter field and the geometry to future outgoing null 
cones, using a nonuniformly spaced set of grid points which follow ingoing null geodesies. 

We find critical behavior at the boundary between black hole formation and dispersion 
for values of the coupling constant 77 in the range of 0.18-100. The critical solution is DSS, 
with the echoing period A strongly depending on rj: As rj tends to 0.18 from above, A rises 
sharply. Moreover we observe small deviations from exact DSS at this smallest 77 value. This 
leads us to conjecture that DSS ceases to be a critical solution for still smaller values of the 
coupling constant. 

The organization of this paper is as follows: In section II we review the basic properties 
of the SU(2) cr-model in spherical symmetry and discuss the system of field equations. We 
present our main physical results in section III, and end the main body of the paper with 
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some conclusions in section IV. In appendix A we discuss our numerical methods, which are 
based on previous work of Goldwirth and Piran (Refs. [10,11]), Garfinkle (Ref. [12]), and 
Gomez and Winicour (Refs. [13-16]). Finally, in appendix B we discuss the convergence of 
our numerical evolutions to the continuum limit as the grid resolution is increased, includ- 
ing both uniform convergence of accuracy diagnostics within a single evolution, and also 
convergence of the numerically computed critical parameter p* itself. 

Conventions are chosen as follows: spacetime indices are Greek letters, SU(2) indices are 
uppercase Latin letters, the spacetime signature is (—,+,+, +), the Ricci tensor is defined 
as R^y = R^xu^ with the sign convention of Ref. [17], and the speed of light is set to unity, 
c=l. 



II. THE SU(2) (7-MODEL IN SPHERICAL SYMMETRY 

Nonlinear cr-models are special cases of harmonic maps from a spacetime (M, g^ u ) into 
some target manifold (N, Gab) (see, e.g., Ref. [2]). Harmonic maps X A (x^) are defined as 
the extrema of the simple geometric action 

S = -f / d?xy/\i\ g^d^X A d v X B G AB (X) . (2.1) 

If the spacetime metric is dynamically coupled to the matter fields X A , then (2.1) must 
be supplemented by the Einstein-Hilbert action. 

Variation of the total action with respect to the a field X A and the metric yields the 
coupled Einstein-a field equations. The stress-energy tensor resulting from (2.1) obeys the 
weak, strong and dominant energy conditions (Ref. [18]). The coupling constant f 2 and the 
gravitational constant G enter the equations only in the dimensionless product 77 = AnGf 2 , 
thereby defining a one-parameter family of distinct gravitating matter models. The field 
equations are scale invariant. 

For the SU(2) u-model, the target manifold is taken as S 3 with Gab the "round" metric 
of constant curvature. Note that the coupling r\ may be interpreted as the inverse of the 
scalar curvature of the target manifold. In the limit i] — > 00 our model thus corresponds to 
the cr-model with 3-dimensional flat target manifold. (This is also easily checked by rescaling 
the field — * (p/^/V an d performing the limit rj — > 00 in Eqs 2.10 and 2.12 - 2.14.) We 
restrict ourselves to spherically symmetric harmonic maps coupled to gravity, which implies 
that the base space (spacetime) must share this symmetry. 

We introduce a Bondi coordinate system {u, r, 9, if} on spacetime based upon outgoing 
null hypersurfaces u = constant, with the line element 

ds * = _ e 2P(u,r) du (Y-^llldu + 2dr^j + r 2 (d6 2 + sin 2 6d V 2 ), (2.2) 

and assume that spacetime admits a regular center r = of spherical symmetry. This 
requires the metric functions near the origin to behave at fixed retarded time uq like 

(3(u ,r) = O(r 2 ), (2.3) 
V(uo,r) = r + 0(r 3 ). (2.4) 
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where the guage has been fixed such that the family of outgoing null cones emanating from 
the center is parametrized by the proper time u at the center. Radial ingoing null geodesies 
are obtained by integrating the equation 

In spherical symmetry the null expansions 0± of inward and outward directed null rays 
emanating from r = constant surfaces, can be defined as Q± = 2(£±r)/r, where £± is the 
Lie-derivative along the null directions l + = e~ 2/3 d r and /_ = 2d u — (V/r)d r . Thus we have 

2 „„ 2 V 

e + = -e~^, 9_ = — -. (2.6) 

r r r 

Whenever + vanishes on some 2-sphere r = constant, this sphere is marginally outer 
trapped. Since this means diverging (3, the Bondi-like coordinate system (2.2) cannot pen- 
etrate a marginally outer trapped surface - in particular an apparent horizon. 

We introduce polar coordinates (0, 0, $) on the target manifold (S 3 , G) , and write the 
SU(2) line element as 

ds 2 = d<p 2 + sin 2 {dQ 2 + sin 2 d$ 2 ) . (2.7) 

We focus on a particular spherically symmetric harmonic map (a corotational equivariant 
map) obtained via the well-known hedgehog ansatz: 

(f>(a?) = <i>(u,r), 0(O = 0, $(:^) = v?. (2.8) 

With this ansatz two of the three coupled fields are determined and only one field <p{u, r) 
enters the equations. Regularity at the origin forces the cr-field to vanish at r = 0, so the 
origin is always mapped to one of the poles of S 3 , defined by the choice of coordinates (2.7). 
As represents the "areal coordinate" of the polar coordinate system (2.7) on the target 
manifold, its regularity behavior near the origin is the same as that of the areal coordinate 
r: 

0K,r)=O(r). (2.9) 
The matter field equations are then reduced to the single nonlinear wave equation 

n<A = ^, (2.io) 

where □ is the wave operator g^ v V ^Vj,: 

D = ^ ((^T + ) d r - - + jd^ . (2.11) 

The nontrivial Einstein equations split up into the hypersurface equations (the {rr} and 
{ur} — (Vy2r){rr} components of = kT^ u ) 
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V = e 2/3 (l-2??sm(0) 2 ) 
the subsidiary equation (r 2 ({uu} — (V/r){ur})) 



V - 2Vf3 = 2r, 



and the redundant equation ({09}) 



(r0) 2 — — (r0') (r<p) 



(2.12a) 
(2.12b) 



(2.13) 



V{rf3" - f3') + rP'V + \rV" - 2r 2 d ur [3 
= rrr<f>'{-V<f>' + 2rd u <P) . 



(2.14) 



The combination of the hypersurface equations (2.12) and the matter field equation (2.10) 
suffices to evolve all the dynamical fields V/r, (3, and 0. Assuming these equations to 
be satisfied, the redundant equation (2.14) then holds identically, and if the subsidiary 
equation (2.13) is satisfied on some r = constant surface (this is assured for r = by the 
regularity conditions there), then it too must hold everywhere. 

In view of this, we construct initial data on a u = constant slice by choosing <p as 
free data on the slice, then integrating the hypersurface equations (2.12) to obtain the 
metric coefficients V/r and j3 on the slice. To evolve this data to future u = constant 
slices, we simultaneously integrate the hypersurface equations (2.12) and the matter field 
equation (2.10). Throughout the initial data construction and the evolution, we use the 
subsidiary equation (2.13) and the redundant equation (2.14) solely to check the accuracy 
of our numerical computations. We discuss our numerical treatment of all these equations 
in appendix A. 

In our coordinates, the Misner-Sharp mass function (Refs. [19-21]) can be written directly 
in terms of the metric, 



r ( V 

m(u,r) = m MS (u,r) = - 1 e~ w 

2 V r 



(2.15) 



or by using the Einstein equations, rewritten as a radial integral within a single slice, 

r 

m(u,r) = m p (u,r) = / m'(u,f)df, (2.16a) 



where 



m'{u,r) = V (-e- 2 P(<P') 2 + 2 
2 \ r 



sin 



(2.16b) 



Since our coordinates would be singular on an apparent horizon, we have designed our 
numerical evolution scheme to slow down as an apparent horizon is approached, in such a 
manner that the evolution only asymptotes to the apparent horizon (cf. appendix Al). 

In other words, none of our numerically-computed slices ever actually contain an apparent 
horizon. Thus strictly speaking we can never measure a black hole mass, but only estimate 
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what the mass will be when (if) a black hole eventually forms. To do this, at each numerical 
time step we compute the Misner-Sharp mass function m MS , and look for regions of the 
numerical grid which are almost at the critical density for black hole formation, i.e. where 
2m MS /r is almost 1. More precisely, if anywhere in the grid 2m MS /r exceeds a specified 
threshold 1 , then we estimate that a black hole will form, with a final mass m BH given by the 
mass function m MS at the outermost such grid point. In general this mass estimate changes 
during the evolution; we use the last value before a numerical evolution terminates as our 
overall estimate for the black hole mass. 

It is also of interest to compute the total mass m to tai within the outer grid boundary. 2 
This gives an upper bound for our final estimate m BU . 



III. RESULTS 



For each value of the coupling constant i], we consider a 1-parameter family of initial 



data 



(no, r), such that (say) for small values of p this initial data eventually disperses 



without forming a black hole, while for large values of p it eventually forms a black hole. 
By using a binary search in p, we can find (a numerical approximation to) the critical value 
p = p* which defines the threshold of black hole formation. 

We have studied the Einstein-a-model system in this manner over the range of coupling 
constants 0.18 < rj < 100, using several different initial-data families. Here we present 
results using the Gaussian-like initial data family 



4>(uq, r) = Ar 2 exp 



r 



(3.1) 



with the "amplitude" A as the parameter p (holding a and r constant for a given critical 
search), and also using the "derivative of 4th-power pseudo- Gaussian" family 



(p{u Q ,r) = -AAr 2 



r — r 



exp 



r — r 



a 



(3.2) 



with the "width" a as the parameter p (holding A and r constant for a given critical 
search). 3 All the results reported here used a "position" r = 5 and an initial-slice outer 
boundary of r oute r = 30. Table I shows some near-critical initial data parameters. 

We have also carried out a number of convergence tests of our numerical scheme, both 
for single evolutions and for entire critical searches. We discuss these in appendix B. 



0.995 for all results reported here. 

2 To be precise, we use m p at the outermost grid point, not m MS , since m MS is numerically somewhat 
ill-conditioned in the outer part of the grid, whereas m p is well-conditioned everywhere. 

3 For this latter case the relative sign of p is reversed with respect to black hole formation, i.e. for 
large p the initial data eventually disperses, while for small p it eventually forms a black hole. 
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A. DSS Echoing 



Discrete self-similarity is denned by the existence of a discrete diffeomorphism $a such 
that for some fixed Ael, 

($* A ) n g = e 2nA g VnGN. (3.3) 

In adapted coordinates r = — In and p = -r^-, where u* is a real number which denotes 
the accumulation time of DSS, we have 

Z(r + nA,p) = Z(r,p) VnGN, (3.4) 

where Z denotes (3, V/r, 0, or any combination of these, e.g. 2m/ r. In addition the a-field 
(f> satisfies the stronger condition 

($i /2 ) n 0=(-ir0, (3.5) 

so that fields even in <fi (e.g. /3, V/r, and quantities constructed from them) are actually 
periodic in r with period A/2. As a DSS diagnostic, we typically look for (A/2) periodicity 
in the black hole formation diagnostic max 2m/ r, where the maximum is taken over r within 
each u = constant slice. 

We have clear evidence for the existence of a type II critical collapse with a DSS critical 
solution. Figure 1 shows examples of this for two values of the coupling constant. Since 
max 2m /r periodicity is only a neccesary condition for DSS, we have also explicitly ver- 
ified that the matter field at selected times u coincides with its image under the DSS 
diffeomorphism <3>; figure 2 shows an example of this. 

We find that the self-similarity echoing period A/2 varies strongly with the coupling 
constant 77. Table I gives some numerical data showing this, and figure 3 shows this same 
data graphically. At large 77, A/2 asymptotes to 0.2300±0.0003. As 77 decreases towards the 
lower limit of the data in table I, rj — 0.18, A/2 rises sharply. At the very smallest coupling 
constant 77 = 0.18, but not at 77 = 0.20 or any larger value, the critical solution shows small 
deviations from exact DSS: the periods and shapes of the individual max 2m/ r oscillations 
deviate by 5-10% from the best-fitting DSS prediction. The physical significance of this is 
not yet clear. 

B. Scaling and Universality 

In the presence of DSS, the black hole mass m BH of slightly-supercritical evolutions shows 
a universal scaling law (Refs. [22-25]) 

In m BH = 7 ln(p — p*) + \1> (ln(p — £>*)) + constant (3.6) 

where 7 sets the overall slope of the scaling law, and the function \l/ is periodic with period 
iA/7 in \n(p —p*). For slightly-subcritical evolutions, the maximum (taken over u within 
each evolution) of the 4-Ricci scalar evaluated at the origin, -R max , also shows a similar 
scaling law, but with slope —27 (Ref. [26]). This is also true for supercritical evolutions, 
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with -R max now defined by taking the maximum in u within each evolution only until a 
(null) slice reaches the apparent horizon. (Our actual evolutions terminate slightly before 
the apparent horizon, but -R max doesn't change significantly in this interval.) 

We have investigated these scaling laws using a sequence of supercritical evolutions with 
varying ln(p — p*). We extract 7 by least-squares fitting In linear function of 

ln(p — p*); after subtracting this fit from lnm BH , we are left with the periodic fine structure. 
Figure 4 shows a typical supercritical scaling law and figure 5 its fine structure. 

Since the numerical resolution of our code is limited by the use of IEEE double precision 
floating point numbers, we expect the errors to blow up for p— p* < 1CT 16 which corresponds 
to \n(p — p*) < 35. This can be seen in figures 4 and 5, and also in figure 8(b) (discussed 
in the next section). For p — p* > —10 deviations from the scaling laws are also apparent, 
demarcating the range of validity of linear perturbation theory. 

We find that the mass scaling exponent 7 varies by at most 5% over the range of r\ we 
have studied, asymptoting to 7 = 0.1185 ± 0.0005 at large 77. (The error is estimated from 
the dispersion in 7 values between fits to critical searches with different initial-data families 
and/or finite difference grid resolutions.) 

The periodicity present in the fine structure of the scaling law (figures 4 and 5) can 
be measured directly. Figure 6 shows a comparison of the measured periods with the 
perturbation-theory prediction |A/7 in ln(p — p*). The agreement is excellent. 

Comparing results for different one-parameter families of initial data, we find that the 
critical behavior is universal at all coupling constants m The critical exponent 7 and the 
echoing period A/2 are the same for all critical searches at a given coupling constant, 
regardless of which initial data family is used. For example, table I shows that 7 and A/2 
are the same (to within numerical errors) for even the very different initial data families (3.1) 
and (3.2). 

IV. CONCLUSIONS 

In this paper we have presented a detailed numerical analysis of SU(2) cr-models coupled 
to gravity in spherical symmetry for a wide range of the coupling constant rj. For 0.18 < i] < 
100 we have evidence of universal critical type II collapse behavior. The critical solution is 
DSS. We have observed both the typical mass scaling at the threshold of black hole formation 
of supercritical initial data, and the corresponding scaling of the scalar curvature for both 
sub- and supercritical evolutions. 

Our numerical results are based on an outgoing-null-cone formulation of the Einstein- 
matter equations, specialized to spherical symmetry (our numerical methods are discussed in 
detail in appendix A), we have carried out thorough convergence tests to ensure the validity 
of our results (see appendix B). Notably, we have demonstrated second order uniform-in-r 
convergence of the error diagnostic 5m (measuring finite differencing errors in the Misner- 
Sharp mass function) for even very-nearly-critical spacetimes. We have also demonstrated 
second order convergence for the initial data's critical parameter p*. To our knowledge this 
is the first time the latter has been reported. 

In the limit of large couplings our model corresponds to the cr-model with 3-dimensional 
flat target manifold. This model has already been studied by Liebling [27], where he consid- 
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ered an additional potential. As this potential does not play a role for criticality we should 
observe the same critical solution for large couplings. In fact our results for both the echoing 
period A = 0.4604 and the scaling exponent 7 = 0.1187 are in good agreement with the 
results reported in [27]. 

While we observe at most a small variation of the critical exponent 7 over the range of 
coupling constants studied, the period A of the DSS depends strongly on the value of the 
coupling constant: as rj tends to 0.18 from above the period increases by more than a factor 
of 2 in the narrow range of 0.18 < rj < 0.3. Also, close to the lower limit we observe small 
deviations from exact self-similarity. 

These observations seem to signal a transition region around the value of rj = 0.18. From 
results of the work on the cr-model in flat space (Refs. [5-7]) it is known that there exists a 
critical (threshold) CSS solution. In a recent paper, Bizoh and Wasserman (Ref. [28]) have 
shown numerically that this solution persists when gravity is turned on, at least up to a 
certain value of the coupling constant. Whether or not the CSS solution plays a role at the 
threshold of black hole formation for small couplings is under current investigation. 

ACKNOWLEDGMENTS 

This work has been supported by the Austrian Fonds zur Forderung der wis- 
senschaftlichen Forschung (project P12754-PHY), NSF grant PHY 9510895 to the University 
of Pittsburgh [S.H.], the Fundacion Federico [M.P. and P.C.A.], and G. Rodgers and J. Thorn 
[J.T.]. We thank Piotr Bizoh for stimulating discussions, and for providing us with research 
results in advance of publication. S.H. also thanks Peter Hiibner for stimulating discussions. 
We thank S. L. Liebling for drawing our attention to the fact that our numerical results for 
large couplings match those of his work in [27]. 

APPENDIX A: NUMERICAL METHODS 
1. Overview 

We discretize the coupled Einstein-matter equations using second order finite differencing 
in r within each u = constant slice, and in u along ingoing null geodesies. Our grid points are 
generically distributed non-uniformly within each slice: On the initial slice we place them 
equidistantly in r between the origin and some finite maximum radius r outcr , but thereafter 
they free-fall in towards the origin along ingoing null geodesies. We always maintain a grid 
point at the origin r = 0; when another grid point reaches the origin we drop the point 
previously at the origin from the grid. 

The choice of freely-falling grid points provides some degree of adaptive grid refinement by 
the focusing of geodesies towards regions of strong curvature. Following Garfinkle (Ref. [12]), 
we also gain additional resolution at late times by explicitly refining our grid by a factor of 
two everywhere in the slice, each time we have lost half of the remaining grid points. Again 
following Ref. [12], for some runs we also manually fine-tune the position of the outermost 
grid point on the initial slice (r outer ), so that this grid point will eventually almost hit the 
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strongest-field region of spacetime. This greatly improves the effectiveness of the factor-of- 
two grid refinements, but this method was not required for the results presented here. 4 

By moving our grid points along null geodesies, the physical domain of dependence 
is automatically contained in the numerical domain of dependence, so our time step is 
not restricted by the usual Courant-Friedrichs-Lewy (CFL) stability limit (Refs. [29,30]). 
However, in order to control time resolution we require (following Refs. [10,11]), that 

(V/r)Au<CAr (Al) 

everywhere in the grid, where C is a constant which we typically take to be on the order 
of unity. The time step Am is thus limited such that grid points fall inwards by no more 
than C/2 grid point spacings per time step. Most of our results reported here were obtained 
with C = 1.5. [Note that for a null-cone evolution similar to ours, but with grid points at 
constant r (Ref. [15]), there is a CFL stability limit, which is in fact just (Al) with C = 2.] 

For rO + sufficiently small, a large value of V/r decreases the time step Au as follows: 
From (2.6) it is clear that for small r @ + the function /3 - which is monotonically increasing 
with r - becomes large (it blows up at an apparent horizon). Furthermore, by (2.15) we 
get V/r = e 2/3 (1 — 2m/ r). Outside of the outermost local maximum of 2m/r both e 2/3 and 
1 — 2m/ r are monotonically increasing with r, and thus so is V/r. If the outer boundary of 
the grid is taken sufficiently far out, this is therefore the location of the maximum of V/r, 
and thus of the most stringent slowdown condition. If Au < 10~ 15 (i.e. close to machine 
precision) the evolution is terminated. 

For the remainder of this appendix, we adopt the usual notation where superscripts 
denote "temporal" (u) levels. Figure 7 shows the typical organization of our grid. All 
discretizations in time (u) and space (r) use nonuniform grid spacings to allow for the free 
fall of the gridpoints and the adaptive time stepping (Al). Our numerical scheme uses the 
geometry fields (5, V/r and (V/r)', and the rescaled matter field ip = r<p. 

Assuming that these fields are known at all grid points on the u = u k and u = u k ~ l 
slices, we determine the fields on the u = u k+l slice as follows: 

• For the innermost 3 non-origin grid points in the u = u k+1 slice, we use a Taylor series 
expansion as described in section A 2. 

• We then sweep outwards over the remaining spatial grids of the u = u k+1 slice as 
discussed in section A 3. 



4 By fine-tuning r outcr in this way, we have also observed DSS in the massless scalar (Klein-Gordon) 
field, with up to 5 echoes visible (in the sense of figure 1). This provides a strong additional test 
of our numerical scheme, since the dynamic range of the DSS is much larger in the Klein-Gordon 
case: A/2 rs 1.73 there (as defined by (3.3)), much larger than the values we find for the cr-field. 
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2. Taylor Expansions near the Symmetry Axis 

The coupled Einstein-matter equations and regularity determine the generic behavior of 
■ip near the origin as 

ip(u + Am, r) = dr 2 + c 2 r 3 + c 2 r 2 Au + O (Aw 4 + r 4 ) . (A2) 

Substitution of this series expansion into the hypersurface equations (2.12) yields corre- 
sponding series expansions for the geometry fields (5 and V/r. 

To determine the geometry and matter fields near the origin on the u = u k+l slice, we 
first least-squares fit the functional form (A2) to the numerically computed ip values at the 
5 innermost non-origin points of the u = u k and u = u k ~ 1 time levels (these points are 
marked by large solid circles in figure 7). This determines the coefficients C\ and c 2 . 

For each of the 3 innermost non-origin grid points on the u = u k+1 slice (these points are 
marked by open circles in figure 7), we first integrate the ingoing null geodesic equation (2.5) 
from u = u k to u = u k+1 , as described below. Then, using the coefficients C\ and c 2 , we 
determine ip at this grid point from the series expansion (A2). Finally, we compute (3, V/r, 
and (V/r)' from their corresponding series expansions. 



3. Integration Schemes 

In order to integrate out from the Taylor series region to the outer boundary, our general 
strategy at each grid point is as follows: 

• We first determine the grid point's r coordinate on the u = u k+1 slice by integrating 
the ingoing null geodesic equation (2.5) from u = u k to u = u k+l . 

• We then determine ip at this grid point using a "diamond integral" scheme of Gomez 
and Winicour (Refs. [15,16,14]). 

• We compute the geometry fields by integrating the hypersurface equations one grid 
point outwards on the u = u k+1 slice. 

For the hypersurface equations (2.12) and the geodesic equation (2.5) we use a second 
order iterated Runge-Kutta scheme (adapted from section 5.2.1, equation (5.6), of Ref. [31]). 
For a generic ODE system dy/dx = f (x, y) the scheme is as follows: 

y k p t c 1 d = y k + Axf(x k ,y k ) (A3a) 
yfe+ i = y k + i Ax ( f (a . fcj yk) + f (a** yj+i)) (A3b) 

While this allows straightforward integration of the hypersurface equations (2.12), the 
geodesic equation (2.5) needs special care: The corrector (A3b) requires evaluating the 
right-hand-side function f at the x h+1 time level. For the geodesic equation this requires 
knowing the field V/r on the u = u k+1 slice, which is not yet computed at the time the 
geodesic integration is done. We thus linearly extrapolate the needed V/r value from V/r 
and (V/r)' values one spatial grid point inwards on the same (u = u n+l ) slice. 
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The matter field equation is integrated using a "diamond integral" scheme of Gomez and 
Winicour (Refs. [15,16,14]). The basic idea is to integrate the nonlinear wave equation (2.10) 
over the null parallelogram E spanned by the 4 grid points N, S, W, and E in figure 7. This 
allows the nonlinear wave equation (2.10) to be written as 

ip(N) = ip(W) + ip(E) - iP(S) 

S 

We evaluate the integral numerically by approximating the integrand as constant over the 
null parallelogram S, with a value which is the average of its values at the grid points W 
and E. This gives second order overall accuracy for ip. 



4. Diagnostics 

Within a single evolution, we use several diagnostics to assess the accuracy of our numer- 
ical computations. We numerically check the satisfaction of the subsidiary and redundant 
Einstein equations (2.13) and (2.14). We also compare the two "different" forms of the mass 
function m MS and m p : These are in fact identical by virtue of the Einstein equations, but 
they are computed in very different ways (via (2.15) and (2.16) respectively), and numer- 
ically they will generally differ by a small amount due to finite differencing errors. This 
difference is a useful diagnostic of the code's accuracy. To this end, we define 

5m(u, r) = mMS - mp (A5) 

TTitotal, init 

where m to tai,init = ^ms( m =0, r max ) is the total mass of our initial slice. 8m is then a dimen- 
sionless diagnostic of how well our field variables approximate the Einstein equations; we 
must have \Sm\ <C 1 everywhere in the grid at all times for our results to be trustworthy. 



APPENDIX B: CONVERGENCE TESTS 

We use convergence tests of the type popularized by Choptuik (Refs. [32-34]) both to 
better understand the performance of our numerical algorithms, and to quantitatively assess 
the accuracy of our numerical results. In particular, it is only through such convergence tests 
that we can be confident our conclusions reflect properties of the continuum Einstein-matter 
equations, rather than numerical artifacts. 

As an example of the convergence properties of our computational scheme, we discuss 
a series of near-critical i] = 0.5 evolutions. We begin by considering the effects of varying 
grid resolutions (specified by the number of grid points N) on the critical parameter p*. 
Figure 8(a) shows these effects for the supercritical mass-scaling law. Notice that the dom- 
inant effect is to simply shift each entire critical curve to a slightly different p*[iV]. Table II 
shows these p*[N] values, and their convergence to a continuum limit (which we denote by 
p*[oo]) as the grid resolution is increased. Notice that the ratios of the successive differences 
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(p*[2N\ —p*[N\)/ (p*[4N\ — p*[2N}) are very nearly equal to 4, i.e. the p* values show second 
order convergence to p* [oo] . 

Besides shifting the effective p*, what other effects does varying the grid resolution have 
on the critical behavior? Figure 8(b) shows the same data as figure 8(a), but plotted using 
the usual logarithmic mass-scaling-law coordinates, and with each grid resolution's data 
plotted using that resolution's own p*[N] value. It is clear that the different resolutions all 
yield the same mass scaling law. 

[In order to get the same mass scaling law at different resolutions, it is essential here to 
use each resolution's own p*[N] value, since figure 8(b) shows the mass scaling law continuing 
down to p — p*[N] values some 10 orders of magnitude smaller than the typical p*[N] shifts 
from one resolution to another. Equivalently, if we did not use each resolution's own p*[iV] 
value in figure 8(b), then the mass scaling law would fail to hold below p — p*[N] ~ 10 -5 
(the typical p*[N] shifts seen in figure 8(a)), whereas by using each resolution's own p*[N] 
value, it actually continues down to p — p*[N] ~ 1CT 15 .] 

We now consider convergence behavior within a single evolution, or more precisely be- 
tween the 3 evolutions whose max 2m/ r time developments are shown in figure 9(a): 

(1) The first evolution uses 8000 grid points, with p = p*[8000] + 10~ 12 , so this evolution is 
just slightly supercritical, by about 1 part in 10 10 . This can be seen in the max 2m/ r 
plot: max 2m/ r first oscillates a number of times, then eventually rises to 1. 

(2) The second evolution uses 16 000 grid points, with the same p as evolution (1). Due to 
the shift in the effective p* with N, this evolution is now subcritical, in fact subcritical 
by a relatively large amount: max 2m/ r oscillates only about half as many times as 
in evolution (1), then eventually decays to zero. 

(3) The third evolution also uses 16 000 grid points, but this time p is adjusted to com- 
pensate for the shift in the effective p* with N: we take p = p*[16 000] + 10~ 12 . By 
construction, this evolution is supercritical again, by the same amount as evolution (1); 
in fact its max 2m /r plot is almost identical to that of evolution (1). 

We use Sm as a diagnostic of our code's numerical accuracy for these evolutions. Fig- 
ure 9(b) shows the convergence of Sm to zero for evolutions (1) and (2). These evolutions 
eventually yield very different spacetimes (one forming a black hole, the other not), but 
here we consider u = constant slices at an early enough time, u = 13.08 (shown by the left 
vertical dashed line in figure 9(a)) that the evolutions have not drifted very far apart yet. 
From figure 9(b) it is clear that 5m is almost precisely a factor of 4 smaller at the higher 
resolution than at the lower one, i.e. 5m shows second order convergence to zero, as expected 
from the construction of our finite differencing schemes. Notice also that this convergence 
is uniform, which is a considerably stronger numerical-fidelity requirement than requiring 
only pointwise or gridwise-norm convergence. 

Now consider a convergence test between evolutions (1) and (3). Because evolution (3) 
adjusts p to compensate for the shift in the effective p* with N, these two evolutions have 
very similar behavior, so we can consider much later u = constant slices and still obtain 
good convergence. For example, figure 9(c) shows the convergence of 5m at the relatively 
late time u = 18.59 (shown by the right vertical dashed line in figure 9(a)). The convergence 
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is (again) very accurately second order. In other words, once we compensate for the shift 
in the effective p* with N, we have excellent - and uniformly pointwise - convergence even 
for evolutions that are very close to critical (p — p* [N] here is about 5 orders of magnitude 
smaller than the p*[N] shifts between the two resolutions), and hence very sensitive to small 
perturbations. 
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TABLES 



Initial Data Family (3.1) Initial Data Family (3.2) 



V 


Parameter is A — 

A* A/2 
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A 




- Parameter 
a* 


is a 

A/2 




7 


0.18 


0.019 523 015 


0.5522 


0. 


1063 


0. 


,003 


1 


.083153 54 


0.5478 


0. 


1028 


0.2 


0.018 942 512 


0.4367 


0. 


1091 





.002 


0. 


.615 31749 


0.4327 


0. 


1150 


0.225 


0.018 241056 


0.3464 


0. 


1207 





.002 


0. 


.651519 42 


0.3472 


0. 


1169 


0.25 


0.017 578 042 


0.3043 


0. 


1173 





.002 


0. 


.688 851 73 


0.3046 


0. 


1173 


0.3 


0.016 392 639 


0.2668 


0. 


1152 


0. 


.002 


0. 


.766 003 44 


0.2675 


0. 


.1146 


0.4 


0.014 534866 


0.2452 


0. 


1132 





.002 


0. 


.929 746 89 


0.2445 


0. 


1139 


0.5 


0.013167 548 


0.2386 


0. 


1152 


0. 


.0015 


0. 


.707 335 37 


0.2386 


0. 


1130 


1 


0.009 528 9751 


0.2314 


0. 


1163 





.0015 


1 


.210 138 07 


0.2313 


0. 


1155 


2 


0.006 809 778 3 


0.2295 


0. 


1179 





.0010 


1 


.064 744 72 


0.2305 


0. 


1167 


5 


0.004 333 205 6 


0.2304 


0. 


1183 
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0. 


.734 344 76 


0.2308 


0. 


1178 


10 


0.003 070 144 2 


0.2293 


0. 


.1186 





.0005 


1 


.318 800 46 


0.2312 


0. 


1182 


100 


0.000 972 589 54 


0.2302 


0. 


1187 


0. 


,0001 


0. 


.631472 58 


0.2311 


0. 


1182 



TABLE I. This table shows two families of near-critical initial data parameters for various 
coupling constants 77. For the Gaussian-like initial data family (3.1), we use the "amplitude" A 
as the parameter p (at a fixed "width" a = 1), with a numerical grid of 16 000 grid points. For 
the "derivative of 4th-power pseudo-Gaussian" initial data family (3.2), we use the "width" a as 
the parameter p (with different "amplitudes" A for different coupling constants), with 8000 grid 
points. For each coupling constant and each family, the table also shows the max 2m/ r echoing 
period A/2 of the near-critical evolution, and the mass-scaling-law critical exponent 7 determined 
for the entire critical search. 



N 


p*[N] 


1000 


0.013156 008 


2000 


0.013164 618 


4000 


0.013166 841 


8000 


0.013167405 


16 000 


0.013167548 



5p* ratio 

) 8.61 xlO- 6 

) 2.22x10- 3.87 

) 5.64 xlO- 7 

> 1.42 xlO" 7 ; 6 * ( 



TABLE II. This table shows the convergence of p* with the finite difference grid resolution N, 
for the Gaussian-like data plotted in figure 8. These evolutions used the same initial data param- 
eters as given in table I. The first two columns give the p* values for the various resolutions N. 
The third column gives the differences 5p* = p* [2N] — p* [N] between consecutive p* values as the 
resolution is doubled, and the last column gives the ratios of consecutive differences. The values 
in the last column are very nearly equal to 4, showing second order convergence of p*. 
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FIGURES 



Part (a) 



Part (b) 



1.0 



1.0 




0.0 1 1 1 1 u 0.0 1 1 1 1 u 

-3 -2-10 1 -3 -2-10 1 

— ln(u*— u) — ln(u*— u) 

FIG. 1. This figure shows DSS echoing behavior in the black hole formation diagnostic 
max 1m I r in near-critical (in this case slightly supercritical) evolutions for coupling constants 
rj = 0.5 (part (a)) and 0.18 (part (b)). Notice the much longer period A/2 of the echoes at 
r\ = 0.18. Although it's not apparent to the eye at the scale of this figure, the 0.18 echoes aren't 
exactly identical: they vary in period and shape by 5-10%. 
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FIG. 2. Snapshots of a near-critical evolution of the SU(2) c-field as a function of lnr for 
r\ = 1.0. The frames are evenly spaced in r = — In u u ^ u (this increases towards the accumulation 
time u*). t increases to the right within each row, then downwards between successive rows. 
Observe that eft is the same in frames in the same column but 2 rows apart; this indicates that <fi 
is periodic in r with period A = 0.46. Also notice that </3 is negated between frames in the same 
column of adjacent rows, i.e. it satisfies the half-period self-similarity condition (3.5). 
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FIG. 3. This figure shows the variation of the near-critical max 2m/r echoing period A/2 with 
the coupling constant rj, for the Gaussian-like initial data family given in table I. Notice the rapid 
rise in A/2 at small rj. The error bar for the 77 = 0.18 point is estimated from the dispersion in A/2 
values when fitting different subsets of echoes in figure 1(b); for larger values of r/ this dispersion 
is negligible. 



1 



• lnrrtBH 

llli? m ax 



X 

B 

5 ^ 



ln(p-p*) 

FIG. 4. Supercritical scaling of the black hole mass m BH , and of -R m ax> the maximum (over 
retarded time u within a single evolution) of the scalar curvature at the origin. The slopes of the 
masses and i? max are +7 and —27 respectively. The scaling fine-structure is clearly visible for 
i? max . Its period is found to be 2.099 which is very close to the value ^A/7 = 2.097 predicted 
by perturbation theory and computed from known values of the critical exponents. This series of 
evolutions was done for 77 = 0.5 using a resolution of 2000 gridpoints. 
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ln(p - p*) 

FIG. 5. This figure shows the fine-scale structure in m BH (shown in figure 4) after subtracting 
a linear fit. For these evolutions we disabled the "2m MS /r > 0.995 detected for N time steps" 
stopping criterion in our code (cf. section A 4), running each evolution until An < 10~ 15 ; in this 
case our final slices' outer grid boundaries almost touched the apparent horizon, so m BH and m to tai 
were essentially identical (within < 1CT 10 of each other). 




FIG. 6. This figure compares the quantity 5A/7, as computed from the echoing in max 2m jr 
and the mass-scaling law, to the period of the oscillations present in the fine-structure of the 
mass-scaling law, which is predicted by perturbation theory to be 5A/7. This has been carried out 
for coupling constants ranging from 0.25 up to 5. All evolutions were done with 2000 gridpoints 
radial resolution. 
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FIG. 7. This figure shows our finite differencing grid. Individual grid points are labelled with 
integers to 6, and their ingoing-null-geodesic trajectories are shown as dotted lines. Grid points 
at the origin are marked with small points. Grid points used in the least-squares fitting procedure 
(cf. section A 2) are marked with large solid circles, while the grid points where the field variables 
are calculated from the Taylor series are marked with large open circles. Grid points where ip 
is already known in the diamond-integral scheme (cf. section A3) are marked with solid squares, 
while the grid point where ip is computed in this scheme is marked with an open square. 
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FIG. 8. This figure shows the convergence of the supercritical mass-scaling law with increasing 
finite difference grid resolution, for r] = 0.5 evolutions. Part (a) shows the mass scaling behavior for 
5 different grid resolutions, plotted on linear scales in both p (here the "amplitude" A) and m BH . 
Notice how the main effect of changes in the grid resolution is to simply shift the entire critical 
curve to a slightly different p*. The actual p* values are given in table II, and all these evolutions 
used the same initial data parameters as given in table I. Part (b) shows the same data plotted 
on logarithmic scales, with each resolution's p values being taken relative to that resolution's own 
p* value. Notice that all the resolutions satisfy the same scaling law, even down to p — p*[N] far 
smaller than the resolution shifts shown in part (a); this is discussed further in the text. 
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FIG. 9. This figure shows the convergence of 5m to zero with increasing grid resolution, for 
near-critical r\ = 0.5 evolutions. Part (a) shows the time development of max 2m/r for each of 
3 evolutions described in the text. Part (b) shows the convergence of 5m between evolutions (1) 
and (2), at a relatively early time. Part (c) shows the convergence of 5m between evolutions (1) 
and (3), at a relatively late time. (Note that in all cases, the marked points are spaced for ease of 
reading, and represent only a small subset of the time steps or spatial grid points.) 
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